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Preface 


This  report  is  the  result  of  a  twelve  week  study  on 
the  feasibility  of  using  the  method  of  weighted  residuals 
to  determine  approximations  to  the  discrete  Green's  function 
or  an  analog  to  it.  The  study  was  sponsored  by  Mr.  Nick 
Pagano,  AFWAL/MLBM.  Included  in  this  report  are  derivations 
of  the  methods  of  Galerkin,  collocation,  and  finite  differ¬ 
ences,  for  the  one  and  two-dimensional  Poisson's  equation. 

The  analytical  solutions  for  several  inhomogeneity  terms  are 
also  presented.  The  results  are  given  in  both  tabular  and 
graphical  form  for  clarity  and  ease  of  reference.  Where  the 
results  showed  a  significant  trend  or  deviated  from  expected 
values,  I  have  attempted  to  provide  an  explanation.  All  of 
the  primary  goals  of  the  study  were  met.  During  the  course 
of  the  twelve  weeks  I  learned  much  about  the  theory  of  Green's 
functions  and  methods  of  numerical  analysis. 

I  wish  to  acknowledge  Dr.  Bernard  Kaplan  for  his  sup¬ 
port  and  direction.  His  valuable  suggestions  aided  me  to 
surmount  difficulties  which  otherwise  might  have  hindered 
the  completion  of  this  study. 

I  must  also  thank  Dr.  Hengehold,  Lt  Col  Bailey,  and 
Mad  Cook  for  their  assistance. 


Randolph  E.  Clapp 
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\  Abstract 


The  purpose  of  this  report  was  to  determine  the  feas¬ 
ibility  of  using  the  method  of  weighted  residuals  to  obtain 
approximations  to  the  discrete  Green's  function  or  an  ana¬ 
log  to  it.  The  methods  of  Galerkin,  collocation,  and  finite 
differences  were  programmed  on  a  CDC  7600  computer  in  For¬ 
tran  IV,  The  resulting  program  was  used  to  generate  the 
approximate  functions  for  the  one  and  two-dimensional  Pois¬ 
son's  equation.  The  two-dimensional  case  was  restricted  to 
the  methods  of  Galerkin  and  finite  differences  on  a  rect¬ 
angular  body.  The  approximate  Green's  functions  and  analogs 
were  applied  to  a  series  of  inhomogeneity  terms  to  obtain 
the  approximate  solutions.  The  results  were  compared  to  the 
analytical  values  at  points  of  interest.  The  average  per¬ 
cent  error  of  the  approximate  solutions  is  reported  for  each 
case  as  the  number  of  interior  nodes  of  the  mesh  was  increas 
ed.  The  areas  of  consideration  were:  the  rate  of  conver¬ 
gence  of  the  approximate  solutions  toward  the  analytical 
solution,  the  amount  of  computer-time  required  to  execute 
tjhe  methods,  and  the  accuracy  of  the  approximate  solutions. 
The  results  of  this  study  indicate  that  the  Green's  func¬ 
tions  and  analogs  obtained  are  valid  approximations  to  the 
discrete  Green's  function  itself,  with  the  restrictions  that 
additional  calculations  may  be  required  in  the  case  of  the 
Galerkin  approximations  and  excessive  computer-time  may  oc¬ 
cur  for  high-accuracy  approximations.  The  finite  difference 
approximations  were  determined  to  be  the  best  method  to  use. 


I.  Introduction 


Background 

In  many  areas  of  scientific  research  one  of  the  major 
steps  in  the  theoretical  modeling  of  experimental  data  is 
the  solution  of  one  or  more  linear  differential  equations. 
Often  these  are  complex,  containing  many  derivatives  and 
source  terms.  Such  equations  are  difficult  to  solve  ana¬ 
lytically.  In  addition,  it  is  sometimes  necessary  to  solve 
large  sets  of  similar  equations  to  determine  trends  and  re¬ 
gions  of  validity  for  theoretical  predictions.  It  is  in 
these  areas  where  computer-aided  approximations  to  the  dif¬ 
ferential  equations  are  useful,  for  they  can  yield  near-cor¬ 
rect  answers  with  only  a  small  amount  of  computer-time. 

Perhaps  the  most  versatile  of  these  approximations  is 
the  Green's  function.  It  allows  the  reduction  of  a  set  of 
linear  differential  equations  to  a  set  of  algebraic  equa¬ 
tions  involving  integrals.  The  solution  set  is  then  found 
by  solving  the  resulting  matrix  equation.  The  exact  ana¬ 
lytical  expression  for  the  Green's  function  is  often  very 

difficult  to  determine. 

% 

In  contrast,  the  discrete  Green's  function,  applicable 
only  at  specific  points,  can  be  approximated  with  much  less 
effort  by  using  finite  difference  methods.  The  solution  set 
can  then  be  obtained  to  the  desired  accuracy  at  these  dis¬ 
crete  points  on  the  body  under  examination. 


Purpose 

The  purpose  of  this  study  is  to  determine  the  feasi- 


bility  of  using  the  method  of  weighted  residuals  (MWR)  to 
obtain  the  discrete  Green's  function  or  a  function  analo¬ 
gous  to  it.  The  methods  of  Galerkin  and  collocation  will 
be  studied  in  comparison  with  the  method  of  finite  differ¬ 
ences.  Both  one  and  two-dimensional  problems  will  be  exam¬ 
ined. 

Plan  of  Attack 

The  approach  to  this  study  will  be  as  follows: 

(1)  To  develop  a  computer  program  which  uses  the 
methods  of  Galerkin,  collocation,  and  finite 
differences  to  obtain  approximations  for  both 
the  discrete  Green's  function,  or  its  analog, 
and  the  solution  for  the  one-dimensional  Pois¬ 
son's  equation.  Homogeneous  Dirichlet  bound¬ 
ary  conditions  will  be  assumed  in  all  cases. 

(2)  To  analyze  the  usefulness  of  the  Green's  func¬ 
tions  or  analogs  obtained  in  the  previous  step 
when  the  inhomogeneity  term  is  varied.  Areas 
of  consideration  will  be:  the  number  of  calcu¬ 
lations  required,  computer  analysis  time,  con¬ 
vergence  rate,  and  overall  solution  accuracy. 

(3)  To  adapt  the  one-dimensional  program  to  handle 
two-dimensional  problems  on  a  rectangular  body. 
Only  the  methods  of  Galerkin  and  finite  differ¬ 
ences  will  be  examined.  Homogeneous  Dirichlet 
boundary  conditions  will  again  be  assumed. 

(4)  To  examine  the  usefulness  of  the  Green's  func¬ 
tions  or  analogs  obtained  as  in  the  above  step. 

«  (5)  To  determine  the  feasibility  and  possible  di¬ 

rections  of  continued  research  into  this  approx¬ 
imation  method. 


Poisson1 s 


Equation  in  One  Dimension 

The  first  problem  to  be  examined  in  this  study  is  the 
one-dimensional  Poisson's  equation.  The  general  form  of 
the  equation  to  be  solved  can  be  expressed  as 

LT(x)  a  g(x)  (1) 


where 


L 

T(x) 

g(x) 


»  The  linear  differential  operator, 


=  The  desired  solution  function 
=  The  inhomogeneity  term 


with  the  associated  homogeneous  Dirichlet  boundary  condi¬ 
tions 


T(0)  =  0  (2-a) 

T(a)  =  0  (2-b) 


Analytical  Solution 

The  general  solution  to  Eq(1)  can  be  determined  by 
direct  integration.  For  the  case  of  a  constant  inhomogen- 
e'ity  the  general  solution  becomes 


T(x) 


C1x  +  0 


2 


(3) 


where  0,  0^,,  and  C2  are  constants.  When  the  boundary  con¬ 
ditions  (Eq(2))  are  applied,  Eq(3)  then  becomes 


The  analytical  solution  for  cases  with  other  inhomogeneity 
terms  can  be  found  in  a  similar  manner. 

Numerical  Approximation 

In  all  of  the  approximation  methods  discussed  in  this 
study  the  body  under  consideration  is  modeled  by  a  mesh. 

The  method  is  then  applied  at  the  N  interior  nodal  points 
of  the  mesh.  This  results  in  a  set  of  N  simultaneous  alge¬ 
braic  equations  which,  when  solved,  yields  the  approximate 
solution  to  the  given  problem.  Depending  upon  which  method 
is  applied  this  approximation  is  valid  for  either  the  entire 
body  or  the  nodal  points  alone. 

The  accuracy  of  the  approximate  solution  depends  on 
the  number  of  nodal  points  in  the  mesh.  A  sample  mesh  for 
the  one-dimensional  problem  with  four  interior  nodes  is 
shown  in  Figure  1.  The  nodes  are  numbered  consecutively 
from  left  to  right.  The  smaller  the  step  size,  h,  the  larg¬ 
er  the  number  of  nodal  points  and  the  more  accurate  the  ap¬ 
proximation. 


Figure  1.  Sample  Mesh  With  Four  Interior  Nodes 


Method  of  Weighted  Residuals.  In  the  general  method 
of  weighted  residuals  (MWR)  the  solution  is  approximated  as 
a  sum  of  N  trial  functions*  fn(x),  weighted  by  undetermined 
factors,  an,  as  (Ref  1:339) 

N 

Vx>  -  lX,fn(x)  (5) 

n=1 


The  trial  functions  are  chosen  to  be  linearly  independent 
members  of  a  complete  set,  each  of  which  satisfy  the  given 
boundary  conditions  (Eq(2))*  According  to  Ozisik  (Ref  Is 
340) ,  the  trial  functions  for  rectangular  coordinates  can 
be  chosen  as  products  of  a  function,  w(x),  with  various 
powers  of  x,  where  w(x)  is  defined  as 


w(x)  -  (x-0)*(x-a) 
«  x(x-a) 


(6) 


It  is  obvious  that  Eq(6)  satisfies  Eq(2)*  For  all  cases  of 
the  one-dimensional  problem  the  trial  functions  will  be  de¬ 
fined  as 

-  fQ(x)  “  xI1*w(x)  (7) 

If  Eq(5)  is  substituted  into  Eq(1)  we  obtain  (Ref  Is 

339) 


L 


r  n 


2“  n,fn(x) 
.  n»1 


g(x) 


(8) 


5 


Multiplying  Eq(8)  by  some  weighting  functions,  6>m, 
integrating  over  the  limits  of  the  body  yields 

a  -a 

u>m(x).Lfn(x)dx  =  I  wm(x).g(x)dx 

0  Jo 


for 


m  =  1 ,2,3*  •  •  • 


Eq(9)  represents  a  set  of  N  simultaneous  equations 
unknowns,  an«  If  we  now  define  the  inner  products 


mn 


<wm,Lfn>  -  (  <J»(x),Lfn(x)dx 

Jo 


®m  =  l  ~m 

f 0 


Jo 


(x).g(x)dx 


Eq(9)  becomes 

N 

2X<‘VLfn>  *  «V«> 

or  in  matrix  form  (Ref  2:6) 

5  5  -  g 


where 


and  then 


(9) 


with  N 


(10-a) 

(10-b) 


(11) 


(12) 


■  denotes  a  square  matrix 
-  denotes  a  column  vector 


,Lf1>  <w1 ,  Lf  2> 
s  <<*>2  ♦  I*f  ,  Lf  2^ 


•  <w1*LfN> 
.  <«2,LfN> 


•  •  •  ^WN*^N^ 


«Aj»S> 

<w2»S> 


«^N»g> 


The  weighting  factors  can  then  be  determined  by 


-  _ 
5  -  I  g 


Finally,  the  approximate  solution  can  be  written  as 


Tn(x)  =  fa 


for  a  method  valid  over  the  entire  body,  or 


T„  -  f« 


for  a  method  valid  only  at  discrete  points,  where 


v  sV*’  ■*- 


W5WW5 


node  V 
Vxnode  2^ 


(20) 


Vxnode  NJ 

f1^xnode  f2^xnode  1^*  •  • 
f1^xnode  2)  *  f2^xnode  2) »  •  * 

f1(xnode  N*»  f2^xnode  *  *  *  *»  fN^*N^ 


•  *  ^jj(xp) 


Green's  Functions  and  Analogs,  Given  a  differential 
equation  of  the  form  of  Eq(1)  and  its  associated  boundary 
conditions  (Eq(2))$  the  analytical  Green's  function  for  the 
problem  can  be  determined.  An  in-depth  analysis  for  the 
case  of  Dirichlet  boundary  conditions  can  be  found  in  Stack 
gold  (Ref  3 :1 -30).  Once  the  Green's  function  has  been  de¬ 
termined  the  solution  to  Eq(1)  with  various  inhomogeneity 
terms  g>j  *62*83*  •••  *6^  can  be  found  by  calculating  the  in¬ 
tegral 

T.(x)  =  (  G(x|x').gi(x)dx'  (22) 

Jo 

where 


G(x|x')  =  The  Green's  function  for  Eq(1)  with 
boundary  conditions  as  in  Eq(2) 


x  =  The  field  point 


x'  =  The  source  point 


For  the  discrete  Green's  function  Eq(22)  takes  the  matrix 
form 


■  SN  “m 


(25) 


where 

?N  i  »  Eq(20)  for  the  inhomogeneity  g^(x) 

GN,jk  “  GN^Xj|Xk^ 

x .  =  The  x  coordinate  of  the  jth  field  node 
J  (d  =  1 *2,3*  • •• »N) 

x.  =  The  x  coordinate  of  the  kth  source  node 
K  (k  =  1,2,3,  • • • »N) 

gN  j(x)  «  A  column  vector  of  order  N  whose  elements 
are  all  equal  to  g^(x) 

*  The  step  size  for  N  interior  nodes 


The  tilde  has  been  placed  over  the  left-hand  side  of  Eq(23) 
to  emphasize  that  the  discrete  Green's  function  solution, 
Hjj,  is  not  necessarily  equal  to  the  weighted  residuals  sol¬ 
ution,  Tjj. 

If  Eq(16)  is  substituted  into  both  Eq(17)  and  Eq(18) 
we  obtain  the  respective  equations 

Tjj(x)  =  f  I"1  g  (24— a) 

and 

TN  -  f  I”1  g  (24-b) 


9 


In  comparing  Eq(23)  and  Eq(24)  it  should  be  noted  that  the 
weighted  residual  approximations  do  not  contain  the  factor 


hN,  This  is  because  the  weighted  residual  solutions  are 
N-term  sumations  whereas  the  discrete  Green's  function  sol¬ 
ution  is  the  approximation  of  an  integral  (Eq(22)).  We  can 
identify  the  analogs  to  the  discrete  Green's  function  in 
Eq(24)  as 

Gn  -  ?  I"1  (25-a) 

for  a  method  valid  over  the  entire  body,  and 
^  — 

Gj,  -  I  I"1  (25-b) 

for  a  method  valid  at  only  discrete  points,  where  the  no¬ 
tation  (t)  is  used  to  denote  the  analogs.  The  functions 
in  Eq(25)  are  termed  analogs  since  the  elements  of  the  ma¬ 
trix  g,  defined  by  Eq(10-b),  in  Eq(24)  are  not  necessarily 
equal  to  g(x). 

Method  of  Galerkin.  In  the  method  of  Galerkin  the 
weighting  functions,  w  ,  are  taken  to  he  equal  to  the  trial 
functions,  f  -  as  defined  in  Eq(7)*  For  these  weighting 
functions  the  inner  products  in  Eq(10)  take  the  form 


„>  ■  r<. 

Jo 


I  =  <f,Lf 
mn  mT 


6m  *  <fm»6>  - 


(x)#Lfn(x)dx 


g(x)dx 


(26-a) 


(26-b) 


The  resulting  form  of  Eq(11)  is  then 


N 

2X<fm-Lfn>  “  <V®> 

n=1 

Since  there  are  no  restrictions  on  the  value  of  x  in  Eq(2?) 
the  method  of  Galerkin  approximations  are  valid  over  the 
entire  body  and  are  given  by  Eq(24-a). 

The  appropriate  form  of  the  analog  to  the  discrete 
Green's  function  is  that  of  Eq(25-a).  It  is  clear  that  this 
must  be  an  analog  since  the  elements  of  the  inhomogeneity 
matrix,  g,  defined  by  Eq(26-b),  are  never  all  equal  to  g(x). 
This  analog  should  be  as  useful  as  the  discrete  Green's 
function  itself  with  the  exception  that  the  inner  products 
to  determine  the  elements  of  the  inhomogeneity  matrix  must 
be  recalculated  for  each  new  inhomogeneity  term. 

Method  of  Collocation.  In  the  method  of  collocation, 
or  point  matching,  the  weighting  functions,  are  defined 
as  (Ref  2:10) 

«„  -  tfx-V  <28> 


where  (x)  is  the  Dirac  delta  function  whose  properties 
are  given  as  (Ref  5:1-30) 


i(x)  « 


x  *  0 
x  =  0 


d(x)dx  = 

0 


0  not  in  (a,b) 
0  in  (a,b) 


(29-a) 

(29-b) 


d(x-x'  )f (x)dx  =  f(x') 


C29-c) 


L 

For  these  weighting  functions  the  inner  products  in  Eq(IO) 
take  the  form 


mn 


<d(x-xm),Lfn>  =  f  d(x-xm)Lfn(x)dx  (30-a) 

Jo 


8m  "  <«x-xm)*8> 


“n<V 


[  i(x-x 

Jo 


m 


)g(x)dx  (30-b) 


8(xm) 


where  (Ref  2:10) 


m  a 

xm  FTT 


(31) 


The  resulting  form  of  Eq(11)  is  then 


N 


Sv-w  -  8<*,) 

n=1 


(32) 


The  value  of  x  in  Eq(32)  is  restricted  to  that  of  the  collo¬ 
cation  points,  xm»  Therefore  the  collocation  approximations 
are  valid  only  at  these  points  and  are  given  by  Eq(24— b). 

The  appropriate  form  of  the  analog  to  the  discrete 
Green's  function  is  that  of  Eq(25-b)*  The  elements  of  the 


inhomogeneity  matrix,  g,  are  the  same  as  the  elements  of 
the  corresponding  matrix  for  the  discrete  Green's  function 
as  in  Eq(23)»  hut  the  approximate  solutions,  TN  and  Tjj, 
are  not  necessarily  the  same  since  is  not  necessarily 
the  same  as  the  analog  form,  f  This  analog  should  be 

as  useful  as  the  discrete  Green's  function  itself  except 
that  the  number  of  terms  necessary  to  achieve  a  given  level 
of  accuracy  may  differ  between  the  two  methods. 

Method  of  Finite  Differences.  According  to  Stackgold 
(Ref  3:1-30)  the  defining  relation  for  the  Green's  function 
for  the  Laplacian  operator  is 

(55) 


with  the  associated  boundary  conditions  for  Eq(1)  and  Eq(2) 


G(0|x'  )  *=  0  (34-a) 

G(a|x' )  =  0  (34-b) 

For  the  discrete  Green's  function  on  a  mesh  with  step  size 
hN  Eq(33)  takes  the  form  (Ref  4:315) 


d2GN(x|x') 


dx 


-1 


hJJ  d(x-x') 


(35) 


The  derivative  term  in  Eq(35)  can  be  replaced  by  a  central 
difference  quotient  as  (Ref  5*6) 


d  Gjj(x|x'  )  Gj^Cx+hjjlxM-SGjjCxIx^+GjjCx-hjjl  x' ) 


Substituting  Eq(36)  into  Eq(35)  and  multiplying  by  hjj  yields 

Gjj(x+hj|  )x'  )-2Gn(x  |x'  )+Gjj(x-hjj  |x'  )  =  h^x-x')  (37) 

Applying  Eq(37)  together  with  the  boundary  conditions  (Eq(3^)) 
to  each  of  the  N  interior  nodes  of  the  mesh  results  in  a 
set  of  N  simultaneous  equations  expressed  in  matrix  form  as 


c  GN  =  hN  IN 


(38) 


where 


C  «  The  coefficient  matrix 

Gn(x1|  x1)  G^(x^|  X2)  •  •  •  G^(x^|  Xjj) 
=  ^N^X2^X/1^  %(x2l  x2)  •  •  •  ^N^x2l 


^ ^ ^  X1 ^  X2^ 


(39) 


«  The  identity  matrix  of  order  N 


The  solution  of  Eq(38)  yields  the  approximate  discrete 
Green's  function  matrix  (Eq(39))*  The  solution  to  the  one¬ 
dimensional  Poisson’s  equation  then  follows  as  in  Eq(23)* 
Computer  Analysis 

The  numerical  approximations  discussed  in  the  previous 
sections  (Eqs(24-a),(24-b),  and  (23))  were  programmed  in  For- 


tran  IV  on  a  CDC  7600  mainframe.  The  resulting  computer 
program  is  on  file  with  the  AFIT  Department  of  Physics, 
Problem  Set.  A  series  of  four  problems  were  ana¬ 
lyzed  as  the  number  of  interior  nodes  of  the  mesh  was  in- 


creased.  The 

areas  of  consideration 

were :  convergence 

rate,  computer-time,  and  the  overall 

The  four  problems  analyzed  were: 

solution  accuracy. 

I) 

IT(x)  «  10 

(40-a) 

II) 

I/T(x)  =  x2 

(40-b) 

III) 

I/P(x)  »  x2  +  1 

(40-c) 

IV) 

IZP(x)  »  x2  +  x  +  1 

(40-d) 

The  Dirichlet  boundary  conditions  for  all  four  problems 
were  taken  as 

T(0)  *  0  (41 -a) 

T(2)  =  0  (41-b) 

The  analysis  was  performed  at  x  ■  0.66  and  x  =  1.33  for  all 

C&S68* 

Exact  Solution.  The  analytical  solutions  to  Eqs(40) 
vtere  found  by  direct  integration  to  be 

I)  T(x)  =  5x2  -  5x  (42-a) 

II)  T(x)  »  ^  “  22T  (42'b) 

III)  T(x)  =  f?  +  f”  "  <*2-e) 

IV)  T(x)  “  %  + 


(42-d) 


The  values  of  the  above  solutions  at  the  points  of  interest 
are  listed  in  Table  I. 


Average  Error.  The  average  percent  error  for  the  one¬ 
dimensional  problem  is  defined  as 


100  I |Tn(0.66)-T(0.66)|  [Tn(1 .33)-T(1 ,33)| 

2  I  T(0.66)  +  T(1 .33) 


(43) 


where 


TN  *  The  approximate  solution  under  consideration 

Convergence  Rate.  The  rate  at  which  each  of  the  three 
approximate  solutions  (Eqs(17) *(18) ,  and  (23))  converged 
toward  the  exact  values  listed  in  Table  I  was  examined  for 
problems  I-IV.  The  average  percent  error  in  the  approximate 
solutions  is  plotted  vs  the  number  of  interior  nodes  in  the 
mesh  as  Figures  2-5.  Appendix  A  contains  the  actual  values 
of  the  approximations  at  the  points  of  interest  for  each  of 
the  four  problems. 


Table  I 

Exact  Solution  Values  at  Points  of  Interest 


Problem 

No. 

Point  1 
x»0.66 

Poinji_2 

x-1.33 

I 

-4.44 

-4. 55 

II 

-0.42798 

-0.62551 

III 

-0.87243 

-1.06996 

IV 

-1.26749 

-1.56379 
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Figure  3.  Average  Error  vs  Node  Number 
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Figure  5.  Average  Error  vs  Node  Number 
for  g(x)  ■  +  x  +  1 


In  each  of  the  four  cases  the  finite  difference  approx' 
imation  yielded  the  best  initial  solution.  In  the  case  of 
problem  I  the  method  yielded  the  exact  values.  Average  er¬ 
rors  below  one  percent  for  eight  or  more  interior  nodes  were 
achieved. 

The  method  of  Galerkin  approximations  converged  at  a 
somewhat  slower  rate  and  did  not  achieve  an  average  error 
below  one  percent  until  20  interior  nodes  were  used.  In 
fact  a  slight  increase  in  the  error  was  seen  at  17  nodes  in 
each  case  examined.  This  was  followed  by  a  decrease  to  be¬ 
low  one  percent  at  20  nodes. 

The  collocation  approximations  yielded  the  worst  ini¬ 
tial  solutions  and  showed  marked  oscillations  in  the  average 
percent  error.  The  oscillations  were  damped  such  that  the 
error  still  dropped  to  around  one  percent  at  the  17  node 
point.  The  slight  increase  and  subsequent  decrease  in  the 
Galerkin  approximations  mentioned  above  may  be  evidence  of 
a  smaller  amplitude  oscillation  there  also. 

The  method  of  weighted  residuals  is  sensitive,  to  the 
choice  of  the  trial  functions.  It  may  be  that  the  oscilla¬ 
tions  in  the  collocation  and  Galerkin  approximations  are  a 
result  of  an  inappropriate  choice  of  trial  functions  for 
the  analysis  of  Poisson's  equation.  Babuska  (Ref  6:241-245) 
states  that  the  trial  functions  as  given  in  Eq(7)  are  numer¬ 
ically  unstable  for  Gaussian  elimination  matrix  inversion. 
That  is,  the  matrix  I  in  the  weighted  residual  solutions  is 
very  nearly  singular  creating  significant  errors  in  the  cal- 


culated  inverse.  Although  Gaussian  elimination  was  not  used 
in  the  direct  matrix  inversion  routine  it  is  reasonable  to 
expect  that  the  same  problem  may  be  present  in  many  inver¬ 
sion  methods.  A  different  choice  of  trial  functions  may  les¬ 
sen  or  eliminate  the  problem. 

Computer  Analysis  Time.  The  quantity  used  as  a  meas¬ 
ure  of  the  computer  analysis  time  was  the  actual  execution 
time  of  the  program.  The  time  taken  by  compilation,  input, 
and  output  was  almost  constant  for  all  computer  runs.  The 
execution  time  in  seconds  is  plotted  in  Figures  6-9  vs  the 
number  of  interior  nodes  in  the  mesh. 

For  all  cases  the  execution  time  was  the  same  for  each 
approximation  method  at  eight  or  less  nodes.  Above  eight 
nodes  the  execution  time  steadily  increased.  This  was  to 
be  expected  since  an  increased  number  of  nodes  requires  ad¬ 
ditional  iterations  of  loops  in  the  program. 

The  method  of  finite  differences  required  the  least 
amount  of  computer  time  (less  than  0.2  seconds  in  all  cases) 
with  the  weighted  residual  methods  requiring  considerably 
more.  It  should  be  noted  that  for  problems  III  and  IV  the 
Gralerkin  method  required  less  time  than  it  did  for  problems 
I  and  II.  This  could  not  be  explained.  The  increased  time 
required  by  the  weighted  residual  methods  can  be  explained 
as  the  result  of  the  additional  matrix  multiplication  re¬ 
quired. 

Overall  Solution  Accuracy.  All  three  approximation 
methods  yielded  solutions  with  an  average  percent  error  of 
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Figure  6.  Execution  Time  vs  Node  Number 
for  g(x)  »  10 
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Figure  8.  Execution  Time  vs  Node  Number 
for  g(x)  »  x2  +  1 
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Figure  9*  Execution  Time  vs  Node  Number 
for  g(x)  «  x2  +  x  +  1 


about  one  percent  when  17  or  more  nodes  were  used  (taking 
the  minimum  point  of  the  collocation  oscillations).  The 
convergence  trends  indicate  that  for  more  than  25  nodes 
the  three  methods  should  yield  nearly  the  exact  values  with 
average  errors  less  than  0.1  percent. 

It  should  be  noted  here  that  for  the  collocation  ap¬ 
proximation  the  direct  matrix  inversion  routine  used  (IMSL 
routine  LINV2F;  Ref  7)  issued  a  warning  that  its  accuracy 
test  had  failed.  This  occured  at  five  or  more  interior  nodes 
Some  of  the  oscillation  in  these  approximations  could  be  a 
result  of  this  round-off  error. 

Conclusions 

From  the  analyses  performed  it  seems  that  for  approx¬ 
imations  with  few  terms  the  method  of  finite  differences  iB 
the  best  choice.  Not  only  is  it  easier  to  apply  (less  cal¬ 
culations  required  and  less  execution  time)  but  the  approx¬ 
imate  solutions  are  more  accurate. 

For  approximations  of  more  than  10  terms  the  trade-offs 
between  calculations  required  and  approximation  accuracy  make 
the  Galerkin  and  collocation  methods  about  equal.  The  finite 
difference  method  is  better  than  both. 


III. 


The  second  problem  to  be  examined  in  this  study  is 
the  two-dimensional  Poisson's  equation.  The  general  form 
of  the  equation  to  be  solved  can  be  expressed  as 


L*T* (x,y)  =  g*(x,y) 


(44) 


where 


1  denotes  a  two-dimensional  function 

i 

L  *  The  linear  differential  operator, 


T*(x,y)  *=  The  desired  solution  function 
g* (x,y)  «  The  inhomogeneity  term 

with  the  associated  homogeneous  Dirichlet  boundary  conditions 


T’(0,y) 

*  0 

(45-») 

t' (a,y) 

=  0 

(45-b) 

t'(x,0) 

-  0 

(45-o) 

T* (x,b) 

*S  0 

(45-d) 

This  is  Poisson's  equation  on  a  rectangular  body. 


Analytical  Solution 


The  analytical  solution  to  Eq(44)  with  boundary  con¬ 


ditions  Eq(45)  can  be  found  by  a  Pourier  series  expansion 
as  (Ref  8:41-42) 


For  a  given  inhomogeneity  the  integration  can  easily  be 
performed  and  the  series  solution  can  then  be  programmed 
to  the  desired  number  of  terms  providing  Eq(46)  converges 
Numerical  Approximation 

A  sample  mesh  and  node  numbering  scheme  is  shown  as 
Figure  10  for  the  rectangular  body.  The  x  and  y  step  sizes 
h  and  k  respectively,  control  the  number  of  nodal  points 
and  hence  the  accuracy  of  the  approximation. 

Method  of  Weighted  Residuals.  In  the  two-dimensional 
MWR  the  solution  is  approximated  as  a  sum  of  N  trial  func¬ 
tions,  f^(x,y),  weighted  by  undetermined  factors,  aQ,  as 
(Ref  1:339) 


N 


TN(x*y)  =  XX  fn(x*y> 
n-1 


(47) 


The  trial  functions  each  satisfy  the  given  boundary  condi¬ 
tions  (Eq(45)).  According  to  Ozisik  (Ref  1:340-344),  the 
trial  functions  for  rectangular  coordinates  can  be  chosen 
as  products  of  a  function,  w  (x,y),  with  various  powers  of 
x  and  y,  where  w'(x,y)  is  defined  as 


w'(x,y)  =  (x-0)*(y-0)»(x-a)»(y-b) 


(48) 


(0,b) 


(a,b) 
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Figure  10.  Sample  Mesh  for  a  Rectangular  Body 
Eq(48)  satisfies  the  boundary  conditions  as  expressed  in 
Eq(45)«  For  all  cases  of  the  two-dimensional  problem  the 
trial  functions  will  be  defined  as 


2n-2 
iT”  ’ 


(xy)  ?  •  w  (x,y)  (for  n«1,4,7,  ...)  (49-a) 


2n+2 


x  *  •  w  (x,y)  (for  n=2,5,8,  ...)  (49-b) 


y  5.  w'(x,y) 


(for  n=3,6,9»  ...)  (49-c) 


If  Eq(47)  is  substituted  into  Eq(44)  we  obtain  (Ref 


Multiplying  Eq(50)  by  some  weighting  functions,  and  then 
integrating  over  the  limits  of  the  rectangle  yields 


a 

<4(x,y)  L  fn(x,y)dxdy 
0 


a 

wm(x,y)  g  (x,y)dxdy 


0 


(51) 


for 


m  —  *1,2,3$  •  •  •  , N 

This  again  represents  a  set  of  N  simultaneous  equations 
with  N  unknowns,  a  .  The  matrix  form  of  Eq(51)  is  the  same 
as  that  for  the  one-dimensional  analysis 


—i  _ i  _i 

I  a  -  g 


(52) 


where 


—  i 
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«  f  1 

i  •  i 

i 

<H,,L  f^> 

,L  fp>  .  .  . 

<4V 

i  i  « 

i  i  i 

f 

<w2,L  £j> 

<w2»L  f2>  •  •  • 

<CJ 

2» 

•  •  9 

1  1  t 

i  «  • 

•  • 

1 

*2>  •  •  • 

<(V 

1 

t 


i  i 

N: 

i  i 


-r 

cr 


«  f 


(55) 


(54) 


<<4pg> 

I  I 

<wo,S> 


fY 

<M,v>  =  |  \  /i*|/dxdy 

JoJo 


The  weighting  factors  can  then  be  determined  by 


a  ■  I  g 


Finally,  the  approximate  solution  can  be  written  as 


I  .  v  211  _1 

TN(x,y)  =  f  a 


for  a  method  valid  over  the  entire  body,  or 


—  i  —i  _i 

TN  -  f  « 


for  a  method  valid  at  discrete  points  alone,  where 


si  r  »  i  i  i“i 

?  . fNj 


TN^xnode  I’^node 
TN^xnode  2’^node  7t) 

|  • 

TN^xnode  N’^node 


^(x2,y2)  f2(x2,72)  •  •  •  ^n^x2*^2^ 


^^(Xjn7jg)  ^2^xN*^N^ 


fN<-xN,yn^ 


(62) 


Green's  Functions  and  Analogs,  Given  a  differential 


equation  of  the  form  of  Eq(44)  and  its  associated  boundary 
conditions  (Eq(45)),  the  analytical  Green's  function  for 


the  problem  can  be  determined.  An  analysis  can  be  found  in 
a  previous  thesis  by  Gallof  (Ref  8s42-43).  The  two-dimen¬ 
sional  form  of  Eq(22)  then  can  be  written  as 


fbf&  »  i 

T.[(x,y)  =  1  1  G*(x|x' ;y|y')*g^(x,y)dx'dy'  (63) 

Wo 


where 


G  (x|x' ; y | y ' )  « 


The  Green's  function  for  Eq(44) 
with  boundary  conditions  as  in 

EqW) 


x  and  y  =  The  field  point  coordinates 


x’  and  y'  =  The  source  point  coordinates 


For  the  discrete  Green's  function  Eq(63)  takes  the  matrix 


form 


?N,i  *  GH*%,i(x*y>,VkS 


(64) 


where 


VV.sVsV.V. 


TN,i  88  E<l(61)  for  the  inhomogeneity  g^(x,y) 

GN,jk  B  GN(xo,xk;y^lyk^ 

x.  and  y .  «  The  x  and  y  coordinates  at  the  jth  field 
3  J  node  (j-1,2,3,  ...  ,N) 

xv  and  y.  =  The  x  and  y  coordinates  at  the  kth  source 
K  *  node  (k=1 ,2, 3»  ...  ,N) 

h«  =  The  step  size  for  R  interior  nodes  along 
the  x-axis 

kg  =  The  step  size  for  S  interior  nodes  along 
the  y-axis 

R*S  =  N 

The  tilde  has  again  been  placed  over  the  left-hand  side  of 
Eq(64)  to  emphasize  that,  as  in  the  one-dimensional  case 

“i 

(Eq(23)),  the  discrete  Green's  function  solution,  TN,  is 
not  necessarily  equal  to  the  weighted  residuals  solution, 

v 

The  two-dimensional  forms  of  Eqs(24)  are  the  approxi¬ 
mate  solutions 

T^(x,y)  =  fl'"1  g'  (65-a) 

and 

TjJ  -  f  T"1  g'  (65-b) 

In  comparing  Eq(64)  and  Eq(65)  it  should  be  noted  that  the 
weighted  residual  approximations  do  not  contain  the  factors 


hR  and  kg.  This  is  because  the  weighted  residual  solutions 
are  N-term  summations  whereas  the  discrete  Green's  function 
solution  is  the  approximation  of  an  integral  (Eq(63)).  We 


can  identify  the  discrete  Green's  function  analogs  in  two 
dimensions  as 

_  I  ~  I  —  I  A 

GN  -  ?  I  (66-a) 

for  a  method  valid  over  the  entire  body,  and 

^  =  FT"1  (66-b) 


for  a  method  valid  at  discrete  points  only.  These  func¬ 
tions  are  termed  analogs  since  the  elements  of  the  matrix 
g  in  the  weighted  residuals  solutions  are  not  necessarily 
equal  to  the  inhomogeneity,  g' (x,y). 

i 

Method  of  Galerkin.  The  weighting  functions,  cd>m,  are 

i 

again  taken  to  be  equal  to  the  trial  functions,  fn,  as  de¬ 
fined  in  Eq(49).  The  resulting  form  of  Eq(51)  is 

N  *b  -a  *b  «a 

l  fm(x>y)L'fn(x,y)dxdy  -  l  l  fn(x,y)g' (x,y)dxdy 

“  Wo  Wo 

for 


m  =  1,2,5,  •••  ,N 

There  are  no  restrictions  on  the  value  of  x  or  y  in  Eq(67). 
Therefore  the  Galerkin  approximations  are  valid  over  the 
entire  rectangle  and  are  given  by  Eq(65-a). 


The  appropriate  form  of  the  discrete  Green's  function 
analog  in  two  dimensions  is  that  of  Eq(66-a).  This  analog 
should  be  as  useful  as  the  two-dimensional  discrete  Green's 
function  itself  except  that  the  double  integrals  must  be  re¬ 
calculated  for  each  new  inhomogeneity  matrix. 

Method  of  Collocation.  In  the  two-dimensional  method 
of  collocation  the  weighting  functions,  o>m,  are  taken  as 
(Ref  2:10) 

"m  -  •<(?-?„)  (68) 

where  i(x)  is  the  Dirac  delta  function  whose  properties  are 
given  in  Eq(29).  The  coordinates  Xjjj  and  ym  are  defined  as 

x„  ■  The  x  coordinate  of  the  mth  interior  (69-a) 
m  node 

ym  =  The  y  coordinate  of  the  mth  interior  (69-b) 
m  node 


For  these  weighting  functions 


mn 


D&  I  8 

iCx-x^j)  i(y-ym)l/f^(x,y)dxdy 

0 


(70-a) 


i  i 


L  VvV 


«(x-xm) i(y-ym)g' (x,y)dxdy 

0 


(70-b) 


The  resulting  form  of  Eq(5'l)  is  then 


N 

"  s’<Vym>  (71) 

The  values  of  x  and  y  in  Eq(71)  are  restricted,  to  the  coor¬ 
dinates  of  the  collocation  points.  Therefore  the  two-dim¬ 
ensional  collocation  approximations  are  valid  at  the  collo¬ 
cation  points  alone  and  are  given  by  Eq(65-b). 

The  appropriate  form  of  the  discrete  Green’s  function 
analog  is  that  of  j;q(66-b).  As  in  the  one-dimensional  case 
these  analogs  should  be  as  useful  as  the  discrete  Green's 
function  itself  except  the  convergence  rate  may  differ  be¬ 
tween  the  two  methods. 

Method  of  Finite  Differences.  The  defining  relation 
for  the  Green's  function  for  the  two-dimensional  Laplacian 
operator  is  (Ref  5:1-30) 

+  =  a(x-x')»(y-y)  (72) 

tor  dy* 

With  the  associated  boundary  conditions  for  Eq(44)  and  Eq(4-5) 
g'(0|x' ;y|y' )  =  0  (73-a) 

G* (a|x' ;y| y' )  =  0  (73-b) 

g' (x|x' ;0| y' )  -  0  (73-c) 

G* (x|x* ;b| y' )  -  0 


(73-d) 


For  the  discrete  Green's  function  on  a  mesh  with  step  sizes 


hjj  and  kg  Eq(72)  takes  the  form  (Ref  4: 315) 


^GR(x|x’ ;y|y' )  +  d2GR(x|x' ;y |y' ) 


hR"  ks"  4(x-x* ) i(y-y' ) 


The  derivative  terms  in  Eq(7*0  can  be  replaced  by  central 
difference  quotients  as  (Ref  5*6) 


a2GR(x|x';y|y«) 


(75-a 


=  GN(x+hR|x*  ;yiy,)-2GN(x|x'  ;y|  y'  )+GN(x-hR|x'  ;y|y') 


dcGN(x|x';y|y') 


Gr(x|x*  »y+kg|  y'  )-2Gn(x|x'  ;y  |y*  )+GN(x|x'  ;y-kg|y') 


Substituting  Eqs(75)  into  Eq(7^)  and  multiplying  by  kg 


yields 


kg^G^x+hj^jx'  ;y|y,)-2(hR2+ks2)GN(x|x'  ;y|y') 


+kg2GR(x-hR|  x' ;y|y' )+hR2GR(xJx' ;y+kg| y» ) 


+hR2GR(x|x»  jy-kgly1)  -  t^kgtfx-x* )  4(y-y' )  (76 


‘  .  *  d  k.'  m  \  C •  k  •  «T.  «T 1  k  •  -  •  •  <4  ■  « 
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Applying  Eq(76)  together  with  the  boundary  conditions 
(Eq(73))  at  each  of  the  N  interior  nodes  of  the  mesh  results 

p 

in  a  set  of  N  simultaneous  equations  expressed  in  matrix 


form  as 


C  GN  =  hpkg  IN 


(77) 


where 


C  «  The  coefficient  matrix 


GN(x1|x1;y1|y1)  #  #  .  GN(x1  |xN;y1  |yN) 

GN(x2,x1;y2|y1)  .  .  .  GN(x2|xN;y2 |yN) 


GN(xN|x1;yN|y1) 


(78) 


IN  =  The  identity  matrix  of  order  N 

The  solution  of  Eq(77)  yields  the  approximate  discrete 

Green's  function  matrix  (Eq(78)).  The  solution  to  Poisson's 

equation  in  two  dimensions  then  follows  as  in  Eq(64). 

Computer  Analysis 
♦ 

The  one  dimensional  computer  program  was  adapted  to 
handle  the  two-dimensional  Galerkin  and  finite  difference 
approximations  on  a  rectangular  body.  The  equations  pro¬ 
grammed  were  Eq(65-a)  and  Eq(64).  The  resulting  computer 
program  is  on  file  with  the  APXT  Department  of  Physics. 

Problem  Set.  A  series  of  four  problems  were  analyzed 
as  the  number  of  interior  nodes  of  the  mesh  was  increased. 


CtKnW 
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The  areas  of  consideration  were:  convergence  rate,  computer 
time,  and  the  overall  solution  accuracy.  The  four  problems 
analyzed  were: 


I)  L*T ' (x,y)  =  10 


II)  L*T ' (x,y)  =  x2 


III)  lV  (x,y)  =  x2  +  y2 


IV)  L  T  (x,y)  =  x2  +  y2  +  x 


(79-a) 


(79-b) 


(79-c) 


(79-d) 


The  Dirichlet  boundary  conditions  for  all  four  problems 


were  taken  as 


T  (0,y)  = 
T ' (2,y)  = 
t'(x,0)  = 
t'(x,2)  = 


(80-a) 


(80-b) 


(80-c) 


(80-d) 


The  analysis  was  performed  at  the  following  four  points  for 


all  cases: 


Point  1  (x  =  0.66,  y  =  0.66) 


Point  2  (x  ■  1.33*  y  -  0.66) 
Point  3  (x  =  0.66,  y  *  1.33) 


Point  4  (x  »  1.33»  7  ■  1.33) 


Exact  Solutions.  The  analytical  solutions  to  Eqs(79) 
were  found  by  integrating  Eq(46)  with  the  respective  inho¬ 
mogeneity  terms.  The  analytical  solutions  are  presented 
as  a  single  function  where  the  inhomogeneity  is  taken  as 

g'(x,y)  *  ax2  +  fa2  +  yx  +  «y  +  t  (81) 

where  a,  yt  i,  and  t  are  constants.  The  solution  is 
then 


The  values  of  the  above  solutiont  for  the  problems  in  Eq(79) 
at  the  points  of  interest  are  listed  in  Table  II.  Each  lis¬ 
ted  value  was  obtained  using  100  summation  terms. 

Average  Error.  The  average  percent  error  for  the  two- 
dimensional  problem  is  defined  as 


Table  II 

Exact  Solution  Values  for  the  Two-Dimensional 
Problem  Set  at  Points  of  Interest 


Problem 

No. 

Point  1 

Point  2 

Point  5 

t 

I 

-2.41 584 

-2.41584 

-2.41584 

• 

II 

-0.20168 

-0.56879 

-0.20168 

i 

III 

-0.40536 

-0.57047 

-0.57047 

t 

IV 

-0.60297 

-0.85364 

-0.77008 

2.41584 


0.56879 


0.73759 


1.02075 


where 


100  V'  IVP°int  *>  -  T  (P°int  *>1 

<en>  *  ~~Tl~  /  .  - 1 -  (85) 

T  (point  i) 


Tm  =  The  approximate  solution  under  considera- 
w  tion 


Convergence  Rate.  The  rate  at  which  each  of  the  two 
approximate  solutions  (Eqs(65-a)  and  (64))  converged  toward 
the  exact  values  listed  in  Table  II  was  examined  for  prob- 

i  « 

lems  I  -IV  •  The  average  percent  error  in  the  approximate 
solutions  is  plotted  vs  the  number  of  interior  nodes  in  the 
mesh  as  Figures  11-14.  Appendix  B  contains  the  actual  val¬ 
ues  of  the  approximations  at  the  points  of  interest  for 
each  of  the  four  problems. 

In  each  of  the  four  cases  the  Galerkin  approximations 


lenrend 


yielded  the  best  initial  solution.  Initial  average  errors 
of  between  one  and  two  percent  were  noted.  In  contrast* 
the  initial  average  errors  for  the  finite  difference  approx¬ 
imations  ranged  from  18  to  27  percent. 

The  Galerkin  approximations  yielded  average  errors  be¬ 
low  one  percent  for  16  or  more  interior  nodes.  The  finite 
difference  approximations  never  had  an  average  error  small¬ 
er  than  four  percent.  Both  methods  converged  at  the  same 
rate  for  10  or  more  interior  nodes. 

The  above  trends  are  inconsistent  with  the  one-dimen¬ 
sional  cases  where  the  finite  difference  approximations  were 
consistently  better  than  the  Galerkin  approximations.  The 
general  trends  should  be  the  same  for  both  one  and  two  di¬ 
mensions.  The  two-dimensional  finite  difference  approxima¬ 
tions  seem  to  converge  rapidly  but  toward  a  value  well  above 
the  exact  solutions.  It  is  suspected  that  there  is  a  minor 
error  in  the  programming  of  the  finite  difference  method 
which  could  explain  this  discrepancy. 

No  oscillations  comparable  to  the  one -dimensional  co¬ 
llocation  cases  were  noted  except  for  a  slight  increase  in 
average  error  for  the  Galerkin  in  problem  I  with  22  nodes 
and  for  the  finite  differences  in  problem  II  with  22  nodes. 

Computer  Analysis  Time.  The  quantity  used  as  a  meas¬ 
ure  of  the  computer  analysis  time  was  the  actual  execution 
time  of  the  program.  The  time  taken  by  compilation*  input, 
and  output  was  almost  constant  for  all  computer  runs.  The 
execution  time  in  seconds  is  plotted  in  Figures  15-18  vs 
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Figure  18.  Execution  Time  vs  Node  Number 
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for  g  (xty)  «  x+  y+  x 


the  number  of  interior  nodes  in  the  mesh. 

In  all  cases  the  Galerkin  approximations  took  more 
time  to  execute  than  did  the  finite  difference  approxima¬ 
tions.  This  is  due  to  the  additional  matrix  multiplication 
required  in  the  Galerkin  method. 

The  trends  in  Figures  15-18  indicate  that  for  more 
than  28  interior  nodes  the  Galerkin  approximations  will  re¬ 
quire  more  than  one  second  of  computer  time  while  the  fi¬ 
nite  difference  approximations  will  require  only  about  0.5 
seconds.  This  is  ah  important  trend  when  a  high-accuracy, 
many-term  approximation  is  desired. 

Overall  Solution  Accuracy.  The  convergence  trends 
indicate  that  both  the  Galerkin  and  finite  difference  ap¬ 
proximations  are  slowly  converging  and  will  not  reach  an 
average  error  less  than  0.1  percent  unless  a  large  number 
(>60)  of  interior  nodes  are  used.  The  Galerkin  method  will 
yield  better  accuracy. 

Conclusions 

From  the  analyses  performed  it  seems  that  neither  of 
the  two  approximation  methods  will  yield  very  accurate  an¬ 
swers  when  few  interior  nodes  are  used.  At  least  16  nodes 
are  required  before  the  Galerkin  method  will  yield  average 
errors  less  than  one  percent.  The  finite  difference  approx¬ 
imations  never  reach  that  level  of  accuracy. 

For  a  large  number  of  interior  nodes  (between  50  and 
100)  the  two  methods  are  about  equal.  The  accuracy  limita¬ 
tions  of  the  finite  difference  approximations  trade-off 


against  the  large  execution  time  requirements  of  the  Gal- 
erkin  approximations.  If  the  suspected  programming  error 
in  the  finite  difference  method  is  located  and  corrected 
then,  as  in  the  one-dimensional  case,  the  method  should  be 
the  best  choice  to  use. 


IV.  Conclusions  and  Recommendations 
Conclusions 

Three  points  may  be  made  regarding  the  use  of  the  me¬ 
thod  of  weighted  residuals  to  determine  approximations  to 
the  discrete  Green's  function.  First,  the  Galerkin  and  co¬ 
llocation  methods  both  yield  analogs  which  are  as  useful  as 
the  discrete  Green's  function  itself.  They  allowed  the  sol¬ 
ution  of  Poisson's  equation  in  one  and  two  dimensions  (only 
the  one-dimensional  case  was  verified  for  collocation)  with 
four  different  inhomogeneity  terms.  The  accuracy  of  these 
approximate  solutions  was  limited  to  the  accuracies  of  the 
methods  themselves.  That  is,  the  solutions  were  the  same 
as  those  obtained  if  the  methods  were  applied  to  each  prob¬ 
lem  without  incorporating  Green's  functions  or  analogs.  The 
only  major  limitation  to  the  method  of  Galerkin  is  that  the 
inhomogeneity  matrix  must  be  recalculated  for  each  respect¬ 
ive  inhomogeneity  term. 

The  second  point  is  the  apparent  sensitivity  of  the 
weighted  residual  approximations  to  the  choice  of  the  trial 
functions.  The  one-dimensional  collocation  approximations 
displayed  a  damped  oscillation  not  present  in  the  finite 
difference  approximations.  Although  this  may  be  a  result 
of  the  round-off  error  (indicated  by  the  failure  of  the  di¬ 
rect  matrix  inversion  routine's  accuracy  test)  in  the  ma¬ 
trix  inversion,  the  fact  that  convergence  was  still  fairly 
rapid  tends  to  say  it  is  not.  The  more  likely  cause  is  an 
inappropriate  choice  of  the  trial  functions  for  the  problems 
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examined. 

Finally,  the  method  of  weighted  residuals  requires 
more  computer  time  than  does  the  method  of  finite  differ¬ 
ences.  This  does  not  become  important  until  large  numbers 
of  interior  nodes  are  used.  In  the  two-dimensional  analy¬ 
sis  it  was  seen  that  a  large  number  (  100)  of  interior 
nodes  was  necessary  to  insure  an  average  error  of  less 
than  0.1  percent.  In  these  cases  the  execution  time  anal¬ 
ysis  indicated  a  rapid  growth  in  the  time  required  to  the 
extent  that  it  may  become  prohibitive  in  terms  of  computer 
resources.  A  large  part  of  this  problem  might  be  eliminated 
by  the  streamlining  of  the  existing  computer  program  and 
using  more  efficient  routines  where  applicable. 


i 


Three  areas  from  this  study  warrant  further  work. 

The  first  is  a  study  of  the  two-dimensional  method  of  co¬ 
llocation.  This  should  be  the  initial  task  before  any  new 
work  is  undertaken.  This  will  allow  a  full  comparison  be¬ 
tween  the  current  analysis  and  any  new  additions.  The  me¬ 
thods  of  least-squares  and  subsectional  basing  (Ref  2:11-14) 
might  also  be  examined. 

The  second  area  is  the  analysis  of  different  trial 
functions  to  see  if  they  affect  the  accuracy  of  the  result¬ 
ing  solutions.  A  possible  choice  would  be  the  eigenfunc¬ 
tions  of  the  Laplacian  operator  (Ref  6:241-245).  Particu¬ 
lar  attention  should  be  paid  to  the  oscillation  in  the  co¬ 
llocation  approximations. 
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The  last  area  for  further  work  is  the  aforementioned, 
streamlining  of  the  existing  computer  program.  Included, 
in  this  would  also  be  the  location  and  correction  of  any 
errors  already  present.  The  matrix  inversion  program  is 
the  routine  considered  the  most  time  consuming.  It  can  pro 
bably  be  replaced  by  a  more  efficient  routine  while  still 
maintaining  the  desired  degree  of  accuracy.  If  a  lower 
accuracy  routine  is  incorporated  attention  should  be  paid 
to  any  induced  round-off  errors. 
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Appendix  A 

Numerical  Approximations  for  the  One-Dimensional 

Problem  Set 

This  Appendix  contains  the  values  of  the  Galerkin, 
collocation,  and  finite  difference  approximations  to  Eqs(40) 
The  values  are  the  output  of  the  computer  analysis  and  are 
tabulated  according  to  inhomogeneity  and  the  number  of  int¬ 
erior  nodes.  The  average  percent  error  is  also  listed  for 
each  set  of  values. 

Note:  The  notation  (*)  indicates  that  for  the  given 
number  of  interior  nodes  the  direct  matrix  inversion  program 
used  indicated  that  it a  accuracy  test  had  failed. 

The  notation  (")  indicates  a  repeated  value  from  the 


line  above  it 


Table  III 

Galerkin  Approximations  for  g(x)  =10 


Number 

of 

Interior 

Nodes 

Point  1 

x  -  0.66 

Point  2 

x  =  1.33 

Average 

Error 

(JO 

2 

-3. 2921 8 

-4.27984 

14.81469 

5 

-4.17238 

-4.38043 

3.78079 

8 

-4.33062 

-4.40286 

1 .74825 

11 

-4.35919 

-4.37570 

1.73239 

14 

-4.34807 

-4.39567 

1 .63283 

17 

-4.35842 

-4.35400 

1.98518 

20 

-4.35830 

-4.41338 

1.31 850 

Table  IV 


p 

Galerkin  Approximations  for  g(x)  =  x 


Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66 

Point  2 

x  -  1.33 

Average 

Error 

<*) 

2 

-0.35117 

-0.61 454 

9.85043 

5 

-0.40985 

-0.62125 

2.45861 

8 

-0.42040 

-0.62274 

1.10697 

11 

-0.42200 

-0.62087 

1.06953 

14 

-0.42096 

-0.62235 

1 .07273 

17 

-0.42157 

-0.62019 

1.17412 

20 

-0.42347 

-0.62817 

0.31427 

14 

17 

20 


Table  V 


Galerkin  Approximations  for  g(x)  =  x2  +  1 


Number 

of 

Interior 

Nodes 

Point  1 

x  -  0.66 

Point  2 

x  *  1.55 

Average 

Error 

(#) 

2 

-0.68058 

-1 .04252 

12.28890 

5 

-0.82708 

-1.05229 

3.42479 

8 

-0.85546 

-1.06305 

1.41104 

11 

-0. 85791 

-1.05844 

1.37050 

14 

-0.85577 

-1.06191 

1.35099 

17 

-0.85741 

-1.05559 

1.53233 

20 

-0.85950 

-1.06951 

0.77352 

Table  VI 

Galerkin  Approximations  for  g(x)=x2+x+1 


Number 

of 

Interior 

Nodes 

Point  1 

x  »  0.66 

Point  2 

x  -  1.35 

Average 

Error 

(*) 

2 

-0.99863 

-1.52538 

11.83411 

5 

-1 .20401 

-1.54885 

2.98185 

8 

-1.24093 

-1.55408 

1.35820 

11 

-1.2471 6 

-1.55766 

1.56347 

14 

-1.24413 

-1.55252 

1.28185 

17 

-1.24642 

-1.54370 

1 .47352 

20 

-1.24919 

-1.56343 

0.73341 

56 
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Table  VII 

Collocation  Approximations  for  g(x)  =10 


Number 

of 

Interior 

Nodes 

Point  1 

x  --  0 . 66 

Point  2 

x  =  1.33 

Average 

Error 

(JO 

2 

-2.96296 

-2.96296 

33.33333 

5 

*-4.191 51 

*-4.69738 

5.69104 

8 

*-3.97101 

* -3. 97101 

10.65219 

11 

*-4.32441 

*-4.56448 

2.70079 

14 

*-4.17154 

*-4.17153 

6.14048 

17 

*-4.36832 

*-4.52071 

1.71439 

20 

*-4.25132 

*-4.23779 

4.49742 

Table  VIII 


O 

Collocation  Approximations  for  g(x)  =  x 


Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66 

Point  2 

x  =  1.33 

Average 

Error 

(JO 

2 

-0.32922 

-0.52675 

19.43228 

5 

*-0.41112 

*-0.64238 

3.31822 

8 

*-0.39642 

*-0.59395 

6.20983 

11 

*-0.41998 

*-0.63352 

1.57490 

14 

*-0.40979 

*-0.60732 

3.57911 

17 

*-0.42291 

*-0.63060 

0.99919 

20 

*-0.41511 

*-0.61152 

2.62186 

Table  IX 


Llocation  Approximations  for  g(x)  =  yr  + 

Number 

Point  1 

Point  2 

Average 

of 

Error 

Interior 

Nodes 

>:  -  0.66 

x  *  1.33 

(#) 

2 

-0.62551 

-0.82305 

25.68956 

5 

*-0.85027 

*-1.11211 

4.38594 

8 

*-0.79352 

*-0.99105 

8.20995 

11 

*-0.85242 

*-1.08996 

2.08141 

14 

*-0.82694 

*-1.02447 

4.73287 

17 

*-0.8597^ 

*-1.08267 

1.32123 

20 

*-0.84024 

*-1.03530 

3.46453 

Table  X 

Collocation  Approximations  for  g(x) 


Number 

of 

Interior 

Nodes 


2 


5 


8 


11 


14 


17 


20 


Point  1 

x  =  0.66 

Point  2 

x  >  1.33 

Average 

Error 

<%) 

-0.92181 

-1.21811 

24.68903 

*-1.20847 

*-1.62280 

6.54321 

*-1.15702 

*-1.45332 

7.88995 

*-1.23948 

*-1.59179 

2.00020 

*-1.20381 

*-1.50011 

4.54813 

*-1.24973 

*-1.58158 

1.26941 

*-1.22243 

*-1.51511 

3.33400 

Table  XI 

Finite  Difference  Approximations  for  g(x)  =  10 


Number  Point  1  Point  2  Average 

of  _  _  Error 

Interior  x  -  0.66  x  *  1.33  (#) 

Nodes 


-4.44444  -4. 


(JO 

0.00000 


Table  XII 


Finite  Difference  Approximations  for  g(x)  =  x£ 


Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66 

Point  2 

x  »  1.53 

Average 

Error 

(JO 

2 

-0.39506 

-0.59259 

6.47743 

5 

-0.41976 

-0.61728 

1. 61819 

8 

-0.42433 

-0.62186 

0.71818 

11 

-0.42593 

-0.62436 

0.40336 

14 

-0.42667 

-0.62420 

0.25776 

17 

-0.42707 

-0.62460 

0.17905 

20 

-o. 42731 

-0.62484 

0.13183 

Number 

of 

Interior 

Nodes 

Point  1 

x  -  0.66 

Point  2 

x  =  1.35 

Average 

Error 

(*) 

2 

-0.83951 

-1 .03704 

3.42506 

5 

-0.86420 

-1.061 73 

0.85626 

8 

-0.86877 

-1.066 30 

0.38079 

11 

-0.87037 

-1.06790 

0.21433 

14 

-0.87111 

-1.06864 

0.13734- 

17 

-0.87151 

-1.06904 

0.09572 

20 

-0.87176 

-1.06929 

0.06971 

lifference 

Table  XIV 

2 

Approximations  for  g(x)  *  x  + 

Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66 

Point  2 

x  *  1.33 

Average 

Error 

(JO 

2 

-1.23456 

-1.53086 

2.35191 

5 

-1.25926 

-1.55556 

0.78504 

8 

-1.26383 

-1.5601 3 

0.26140 

11 

-1.26543 

-1.56173 

0.14713 

14 

-1.26617 

-1.56247 

0.09428 

17 

-1.26658 

-1.56287 

0.05585 

20 

-1.26682 

-1.56311 

0.04817 

Numerical  Approximations  for  the  Two-Dimensional 

Problem  Set 

This  Appendix  contains  the  values  of  the  Galerkin  and 
finite  difference  approximations  to  Eqs(79)*  The  values 
are  the  output  of  the  computer  analysis  and  are  tabulated 
according  to  inhomogeneity  and  the  number  of  interior  nodes. 
The  average  percent  error  is  also  listed  for  each  set  of 
values. 

Note:  The  notation  (**)  indicates  that  for  the  given 
number  of  interior  nodes  the  direct  matrix  inversion  program 

=i 

indicated  that  the  matrix  I  was  algorithmically  singular. 


The  resulting  solution  values  were  discarded  as  being  invalid 


Table  XV 

Galerkin  Approximations  for  g  (x.y)  «  10 


-2.458191  -2.4363: 


-2.41038  -2.43 


Bill 


-2.39441  -2.4241 


II 


Table  XVI 

Galerkin  Approximations  for  g  (x,y) 


Number  Point  1  Point  2  Point  3  Point  4  Average 
of  Error 

Interior  x  *  0.5S  x  =  1.55  x  ■  0.66  x  =  1.55  (%) 

Nodes  y  *  0.5S  y  »  0.65  y  *  1.55  y  *  1.33 


PI 


-0.19932  -0.370051  -0.200081  -0.38251|  1.50636 


-0.1991 


-0.19927 


0.200141  -0.36458  1.41233 


-0.37342  -0.20255  -0. 36642]  0.94309 


-0.366541  0.71061 


-0.36653  0.6620 


Table  XVII 

»  o 

Galerkin  Approximations  for  g  (x,y)  *  x 


Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66 
y  =  0.55 

I'oint  2 

x  =  1.33 

.7  =  0.56 

Point  3 

X  a  0.55 
y  =  1.33 

Point  4 

x  =  1.53 
y  =  1.33- 

Average 

Error 

(*) 

4 

-0.39863 

-0.57013 

-0.57013 

-0.76502 

1.25049 

10 

-0.39831 

-0.57813 

-0.57813 

-0.72917 

1.27232 

16 

-0.39755 

-0.57597 

-0.57597 

-0.73285 

0.97297 

22 

-0.39854 

-0.57373 

-0.57372 

-0.73308 

0.73910 

28 

-0.39905 

-0.57364 

-0.57361 

-0.73304 

0.70008 

34 

*  * 

*  * 

*  * 

♦  * 

40 

*  * 

*  * 

*  * 

+  * 

Table  XVIII 


i  p 

Galerkin  Approximations  for  g  (x,y)  =  x  + 


y2  +  x 


Number 

of 

Interior 

Nodes 

Point  1 

x  =  0.66) 
y  =  0.55 

Point  2  Point  3 

x  =  1 .33  x  »  0.66 

y  a  0.66  y  -  1.33 

Point  4 

x  =  1.33 
y  =  1 .33' 

Average 

Error 

(50 

4 

-0.59706 

-0.85163  -0.77124 

-1 .05504 

1.18138 

10 

-0.59666 

-0.86734  -0.77552 

-1.01021 

1.09759 

16 

-0.59570 

-0.86139  -0.77607 

-1.01482 

0.86809 

22 

-0.59694 

-0.85854  -0.77330 

-1.01512 

0.63594 

28 

-0.59761 

-0.85842  -0.77328 

-1.01507 

0.64254 

34 

** 

*  *  *  * 

*  * 

40 

*  * 

*  *  *  * 

*  * 

Table  XIX 

Finite  Difference  Approximations  for 


Number  Point  '1 
of  _ 

Interior  x  =  0.66 
Nodes  y  »  0.66 


oint  2 1  Point  3 
x  *  0.66 


Point  4 


x  =  i.22 

y  =  1.35 


Average 

Error 

(*) 


-2.05128  1-1.62393  1-2.13675 


-2.26123  -2.18765  1-2.27367 


2.29383  1-2.26968  1-2.29771 


2.30374  1-2.29305  1-2. 3054-3 


2.30789  (-2.30225  1-2.30877 


2.30993  (-2.30665  1-2.31050 


2.31116  (-2.30905  1-2.31150 


2.05128 


2.25634 


2.29182 


2.30279 


2.30737 


2.30967 


2.31097 


8.56088 


7.03412 


5.20250 


4.66425 


4.44396 


4.33500 


4.27410 


Table  XX 


Finite  Difference  Approximations  for  g  (x,y) 


Number 

Point  1 

of 

Interior 

x  =  0.66 

Nodes 

y  *  0.66 

0,66 

1.55 


Error 

(#) 


0.14815  hO. 23457  1-0.16049  1-0.29630  |  25.75409 


0.18346  1-0.32431  hO. 18583  1-0.33749  !  9.36034 


0.19016  (-0.34152  1-0.19096  1-0.34612  |  6.14224 


0.19232  1-0.34689  1-0.19269  1-0.32858  |  6.48503 


0.19326  1-0.34912  1-0.19345  1-0. 35025  I  4.65414 


0.19373  1-0.35022  1-0.19385  1-0.35089  I  4.42834 


083  1-0.19408 


0.19401 


vW 


Table  XXII 


i  p  2 

Finite  Difference  Approximations  for  g  (x,y)  *  x  +  y  +  x 


Number 

Point  1 

of 

Interior 

x  »  0.66 

Nodes 

y  *  0.66 

Point  3  Point  4  Average 
_  Error 

x  »  0.65  x  =  1.33  (%) 

y  -  1.55  y  =  i .55 


24.62142 


10.23579 


6.94492 


6.61934 


6.45710 


6.36578 


sV.'AK  a'.Va* * 
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